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Abstract. A method for the discrete particle simulation of of almost rigid, sharply 
edged frictional particles, such as railway ballast is proposed. In difference to 
Molecular Dynamics algorithms, the method does not require knowledge about the 
deformation-force law of the material. Moreover, the method does not suffer from 
numerical instability which is encountered in MD simulations of very stiff particles. 

1 Introduction 

Traditionally, the subgrade of railway tracks is modeled using continuum me- 
chanics. These methods have been proven to yield reliable results in many 
applications and have been developed to standard methods. In certain appli- 
cations, however, continuum models fail in describing the mechanical proper- 
ties of the subgrade. This is the case when the ballast must not be considered 
as a continuous medium, but when the granularity of the material is impor- 
tant. Typical processes which cannot be sufficiently explained by continuum 
models are sedimentation of the ballast due to recurrent load, abrasion of 
the ballast particles which leads to less efficient damping properties and the 
formation of force chains inside the ballast material. 

In the past decade mainly by physicists much work has been done in the 
field of Molecular Dynamics of granular material, i.e., the numerical sim- 
ulation of granular material as many-particle systems. This technique was 
applied to many interesting systems and has contributed to the explana- 
tion of several exciting and technologically important effects, such as mix- 
ing and demixing of granular materials, avalanche statistics on sand heaps, 
milling processes, convection dynamics in shaken granular materials and oth- 
ers. Many examples of such simulations can be found, e.g., in Q. 

The idea of Molecular Dynamics is to simulate the granular material as a 
many particle system and to determine the dynamics of the system by numer- 
ical integration of Newton's equation of motion for each of the N particles: 

ta =F i /m i 

.: :' (i) 

& =jr 1 M l , 

where T"j and (pi are the position and the orientation of the i-th particle of 
mass m,i and moment of inertia Ji while Fi and Mi are the force and the 
torque acting on this grain. In three dimensions Eqs. (TO) establish a set of 
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6N coupled non-linear differential equations which have to be numerically 
integrated. The force F{ consists of gravity and of the interaction force of 
the particle i with other particles j 

Fi = m l9 + Fij , (2) 

3 

where g is the gravitational acceleration. Granular particles interact with 
each other only if they are in mechanical contact. For spheres of radii Ri and 
Rj we write 

F _ / F$ n tJ + F T t if \n -r 3 \<R t + Rj 

F »-\0 else, (3) 

where F N and F T are the components of the force in normal and tangential 
direction with respect to the inter-center vector — rj and n and t are 
the corresponding unit vectors. There exist several models for the interaction 
forces in normal and tangential direction F^ T (ri, rj, fj, rj) (see, e.g. 0) 
which shall not be discussed in detail here. 



2 Molecular Dynamics fails for the simulation of 
railway ballast 

There are many examples where granular systems have been simulated by 
Molecular Dynamics, however, except for few examples, realistic simulations 
have been achieved only for systems where the dynamical behavior of the 
grains dominates the system properties. When the static properties of the 
particle system become important, i.e., when the relative velocities are small 
or zero, the details of the interaction force become essential for understanding 
the system behavior. We are faced with two major problems: 

• The interaction force of contacting particles must be known as a function 
of the particle positions and velocities. In the case of sharply edged grains 
such as railway ballast, this function is unknown. 

• As soon as the realistic simulation of static properties matters for the sys- 
tem behavior, the simulation slows down extremely which implies that 
to achieve affordable computation times one has to make simplifying as- 
sumptions on the particle contact which are not justified from the point 
of view of mechanics and material science. 

For these reasons we believe that the described Molecular Dynamics method 
is in principle unsuitable for the simulation of railway ballast. Below we list 
arguments which support this statement: 
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1. Elastic normal force: Whereas the normal elastic force of contacting 
spheres is given by the Hertz law 

f n,oi _ 2 yygjff 3/2 ( .s 

ij 3 1 — v 2 ■? 

with Y, u, and R c g being the Young modulus, the Poisson ratio, and 
the effective radius, respectively, this force is not known for more com- 
plicatedly shaped particles. For smooth particles (when the local ra- 
dius of the contact area is large as compared to the compression £jj = 
Ri + Rj — \r% — rj\) the function (^) is certainly a good approximation, 
for sharply edged particles, however, this function fails. Some authors 
assume that the normal force is proportional to what they call "overlap" , 
i.e., the volume of the compressed material or (yet more simple) in two 
dimensions the "compression length" (for spheres the value £) which is 
certainly incorrect even for the most simple case of contacting spheres 
and the more for more complicatedly shaped grains. 

2. Dissipative normal force: The dissipative normal force of contacting 
bodies is unknown, in general. For viscoelastically colliding spheres and 
other smooth contacting surfaces it is given by 

F N,diss = A ^9_ p N, e l (5) 

The prefactor A is a complicated function (for details see Hj) which con- 
tains the viscous constants of the material which are unknown in general. 
For particles which are not smooth, such as railway ballast, even the func- 
tional form of the dissipative force is unknown. The frequently applied 
force law F^' dlss oc £ is not justified and fails even for spheres. 

3. Tangential force: Whereas for smooth bodies the normal force can be 
determined from bulk properties of the material, the tangential force is 
determined by the bulk and by surface properties. A natural (phenomcno- 
logical) assumption is 

Fl < i*Fg , (6) 

where fi is the Coulomb friction parameter. Unfortunately, this model is 
not sufficient for static systems: Assume a (non-spherical) particle which 
rests on an inclined plane (angle a). Its normal force is F N = mg cos a 
and there is a corresponding tangential force too. To prevent the particle 
from sliding along the plane one has to assume an additional force which 
mimics static friction, e.g. Jl. This force cannot be derived from material 
or surface properties, hence, it is arbitrary. Choices for this function which 
can be found in literatures are even in disagreement with basic mechanics. 

4. Numerical integration of Newton's law: Besides the fact that the 
numerical integration of the set (|l|) for a relevant number of particles 
N, e.g., N = 1,000 which establish a quite small system of 10 x 10 x 
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10 particles in three dimensions, requires serious numerical effort, there 
are principal problems when applying Molecular Dynamics to a particle 
system: 

• For very rigid particles the gradient of the force is large and, therefore, 
requires a very small integration time step. The more rigid the particle 
material the slower progresses the simulation. Assuming a time step 
of At = 1CP 7 sec which is used in many MD simulations of granular 
material, one needs 10 9 integration steps of the equations of motion 

to achieve a real time of 100 sec only. Hence, a simulation of long 
time behavior using MD seems to be unrealistic. 

• The accuracy of the algorithm cannot be increased too much by 
reducing the time step. Frequently, predictor-corrector algorithms 
are used, which require powers of the time step, (At) 2 = 10~ 14 , 
(At) 3 — 10~ 21 , etc. (again for At = 10~ 7 ). If we ignore the prefactor 
of this powers in the integration scheme to obtain a crude estimate, 
one has to sum numbers of the orders 0(1), O (10 -7 ), O (l0~ 14 ), 
etc. Double precision numbers have a mantissa of 15 digits, therefore, 
one cannot sum numbers which are different by more than 10 15 . This 
limits the minimal value of the integration time step. If one requires 
higher precision one needs real numbers in quadruple precision or 
higher which implies higher memory consumption and (more impor- 
tantly) yet larger computation time since multiplication of two num- 
bers in quadruple precision requires approximately quadruple time as 
a multiplication of double precision numbers. 

We are aware that these numbers are a very crude estimate which 
have been given to illustrate the problem. Of course, more sophis- 
ticated integration schemes suffer less critical from the addition of 
different numbers. In principle, however, the problem persists. 

3 Rigid Body Dynamics 

In Molecular Dynamics simulations the trajectories of the particles are deter- 
mined by numerical integration of Newton's equations of motion. As discussed 
above, this implies the "soft particle assumption", i.e., the particle deform 
under load. Since the particles are very hard but not completely rigid we 
encounter serious numerical problems as described above. 

The Rigid Body Dynamics originates from the opposite idea: The inter- 
action forces are determined from the required behavior of the particles. This 
method is, therefore, suited to simulate perfectly rigid particles without the 
necessity to specify a certain force-deformation law. We consider the example 
of a rigid sphere which rests on a rigid flat surface (see Fig. |Tj). There is a 
point-like contact of the sphere and the surface, hence, there is a contact 
force F N in vertical direction. Moreover gravity acts on the sphere causing a 
force —mg. If one would chose the contact force F N = 0, the sphere would 
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2=0 

Fig. 1. A rigid sphere resting on a rigid flat surface. 

move downwards with an acceleration <?, i.e., it would penetrate the surface. 
This unphysical behavior has to be avoided by the proper choice of the force 
F N due to the following conditions: 

1 . Contact forces have to be chosen in order to avoid penetration of contact- 
ing particles. 

For our example F N > mg follows. On the other hand, a force F N > mg 
would lead to unphysical behavior since the sphere would move upwards. 
This is avoided by the second condition. 

2. A contact force vanishes when the contact breaks. 

A contact is said to break if the normal component of the relative accel- 
eration of the concerned particles, or their normal velocity is larger than 
zero. (The relative velocity is counted positive if the particles separate 
from each other.) 

3. There are no attractive normal forces. 

In our example F N > mg causes the contact to break which implies 
that the contact force vanishes. Therefore, from these conditions follows 
F N = mg. For this choice the total force is zero and the sphere rests 
on the plane. The conditions 1.-3. arc sufficient to describe any particle 
system as long as there are no friction forces. For systems with friction 
we need one more condition: 

4. Friction forces act in parallel with the contact plane, i.e., perpendicular to 
F N . Given the tangential force F* which is necessary to keep two particles 
from sliding. Then the acting tangential force is \F\ — min (\F*\ , J^t-F^l). 
Its sign has to be chosen opposite to the tangential relative acceleration 
(or the tangential relative velocity). 

In agreement with Coulomb's friction law the particles slide only if \F T | > 
\F N \. If this condition is fulfilled the friction force adopts its maximal 
value ±fiF N . 

To perform simulations one has to derive the forces of the particles in nor- 
mal and tangential directions from these four conditions. The corresponding 
algorithm will be explained in sections || and [| 

In our simple example (Fig. |l|) there exists only one contact between the 
sphere and the plane. In more complex situations, e.g. for a resting cube, 
there are contacts areas instead of points. These contacts may be always 
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reduced to point contacts. It will be shown that the described conditions 
are sufficient to determine the forces and torques which act on the particles, 
provided there are not too many contacts in the system. If the number of 
contacts is too large only the total force and the total torque which act on 
a particle may be determined, but not each of the pairwise contact forces. 
For the computation of the particle trajectories, however, the total forces and 
torques are sufficient. We will return to this issue in Sec. |^. 

Rigid Body Dynamics has been intensively studied in the past two decades. 
Descriptions of the algorithm can be found in j^H . The core of the algorithm 
is the numerical computation of the contact forces which is a Linear Com- 



plementarity Problem [ j|J7|. An efficient algorithm for this type of problems 
can be found in ]8||j 10,11]. Rigid Body Dynamics has also been applied to 
granular systems, e.g. ]^,[p|, where frictionless smooth spheres have been 
simulated. Systems of granular particles subject to friction have been stud- 
ied, e.g., in Ipf, 

Due to our understanding the Rigid Body Dynamics is much better suited 
for the simulation of railway ballast for the following reasons: 

• Ballast particles are irregularly shaped and sharply edged. Even if the 
bulk material properties were precisely known, the contact force law is 
unknown due to the complicated shape. 

• Ballast particles are very stiff which implies that the gradient of the 
interaction force is very steep. In this regime the numerical integration 
of Newton's equation is problematic. For the dynamics of the system 
the deformation of single particles is unimportant, i.e., the Rigid Body 
assumption is well justified. 

• Static friction, whose treatment in MD-simulations is problematic too, is 
essential for the dynamics of the system. It is correctly modeled in Rigid 
Body Dynamics. 

• The long time behavior of railway ballast is affected by abrasion and 
fragmentation of particles. Molecular Dynamics of fragmenting particles 
may cause artifacts for several reasons which cannot be discussed here in 
detail (see jlj|). Rigid Body Dynamics is very well suited for this case. 



4 Schedule of Rigid Body Simulations 

The state of the granular system is described by the position and orientation 
of its particles and by the according time derivatives. Contacts between the 
particles may be classified as sticking and sliding contacts. In the due of time 
the contact network is modified by creation and termination of contacts as 
well as by transformation of sticking contacts into sliding ones and vice versa. 
Whenever the contact network is modified the state of the system changes 
qualitatively. The simulation proceeds in discrete time steps. Each of them 
consists of 
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1. Contact detection: all existing contacts are registered. 

2. Treatment of collisions: A collision takes place if two contacting parti- 
cles move with negative normal relative velocity. In this case one cannot 
determine a finite contact force which would avoid penetration of the par- 
ticles since any force, however large it is, would need a short but finite 
path to retard the colliding particles. Hence, mutual penetration would 
be unavoidable. Therefore, we need a special treatment for collisions (see 
Sec.0). 

3. Formulation of the geometry equations: After a collision, in general, 
there is a number of contacts of particles which have a positive normal 
relative velocity, i.e., the particles lose contact. These contacts have to be 
erased from the list of contacts. The normal components of the relative 
velocities at all remaining contacts are zero. We have to establish the 
geometry equations which contain the information about the geometry 
of the system (see Sec. |J). 

4. Computation of the forces: Section ^ deals with the computation of 
the relative accelerations by means of the geometry equations. 

5. Integration of the equations of motion: Finally we have to inte- 
grate the equations of motion for all particles. During this operation it 
may be necessary to update the geometry equations and to repeat the 
computation of the forces according to the integration scheme used. 

5 Mathematical description 

For convenience at first we will restrict to frictionless particles. When the 
mathematical framework has been developed for this simplified case we will 
then introduce friction forces between particles. The rigidity of the particles 
is enforced by means of mathematical motion constraints of the form 

g(q) > o , (7) 

where q is a vector containing the positions (center of mass position and 
orientation) of all particles of the system. The constraint function g shall be 
zero if particles in the system are in contact and larger than zero otherwise. 
For spheres, which is the most simple case, the constraint function reads 

g(q) = \r i -r j \-R i -R j . (8) 

If the spheres would deform each other (|rj — rj\ < Ri + Rj) the function 
g(q) would be negative, if the particles touch each other it would be zero. To 
prevent the deformation we require g{q) to be positive or zero, i.e., it has to 
fulfill the condition (0). For sharp edged particles, such as particles which are 
described by polyhedrons or polygons, the motion constraints are 



g(q) = n (n + x. t - r 3 ) - d , 



(9) 
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where rij is normal of the face of particle j which is in contact with an edge 
of the particle i. The vectors 7% and rj are the center of mass positions of 
the particles, the edge of particle i that is in contact with the particle j is at 
position Ti + Xi. The constant d is the distance of the contacting face of j 
from the center of the particle. The left picture in Fig. ^ shows a sketch of a 
typical contact of two particles. Face-face contacts can be described by two 
face-edge contacts (see Fig. 0). 



Fig. 2. Left: a face-edge contact, ri and rj are the center of mass positions, Xi is 
the coordinate of the contacting edge relative to the center of mass of particle i and 
rij is the normal of the contacting face of particle j. Right: face-face contacts can 
be reduced to two face-edge contacts. The face normals of each of the contacts are 
displayed as well. 

Every motion constraint corresponds to a scalar contact force /. Accord- 
ing to d'Alembert's principle the direction of the contact force is given by 
the spacial derivative of g with respect to all components of the coordinate 
vector q, namely dg/dq. The contact force that acts on a certain particle i 
is fdg/dqi, with qi being the coordinates of particle i. We can formulate the 
equation of motion for the particlesn 







(10) 



9a(q) > . 



(11) 



Mi is the mass matrix of the particle i, which has the form 



Mi = 



(mi 0\ 

TYli 

m, 
V OjJ 



(12) 



1 Particle indices are written in Latin letters, contact indices in Greek letters. 
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where Ji is the moment of inertia tensor. In two dimensional systems Ji is 
only a scalar Ji and there are only 2 entries of mi. Qi finally is the external 
force (and torque) acting on particle i. This is usually gravitation, but other 
external forces can be incorporated at this point as well. 

Although we formulated the motion constraints in the form g(q) > 
to allow separation of particles, contact forces can only act if particles are 
actually in contact. Therefore, constraints which are strictly positive, i.e. 
g(q) > (the particles are separated), can be disregarded. These constraints 
are said to be inactive. The remaining constraints, the active ones, are thus 
satisfied by g{q) = 0. Since the g(q) have to remain non-negative their time 
derivatives have to be non-negative as well. We therefore have 

9a = ^^qk > (13) 
oqk 

dg a .. . d 2 g a . . . 
9a = -^—qk + a » qkQl > . 14 
dq k oq k dqi 

For simplicity of notation we used the Einstein convention, i.e. summation 
over doubly occurring indices k and I is implied. These time derivatives are 
the relative velocity and relative acceleration of the particles at their contact 
points. It is important to note that g a and g a are not the relative velocity or 
acceleration of the centers of mass of the particles but of the points of both 
particles which are actually in contact. It can easily happen that the relative 
velocity or acceleration of the contact points are positive (the particles are 
about to separate) although their centers of mass approach each other. 



We insert the equation of motion (10) into (fL4J) and find 




IkOqi 



The first term on the right hand side describes the action of the external 
forces, the second term describes the action of inertial forces as, e.g., cen- 
trifugal force and Coriolis force, the third term finally describes the action of 
the contact forces. We can rewrite this equation as 

9a = b a + ^2 A <*pU i ( 16 ) 



with A a p and b a abbreviating 

A - d9a fcf-1 dg $ 

A-a/3 — "5 M k — 

dq k dq k 



(17) 

dqk V"" fc """7 ' dqkdqi 



°9a f^r-l^ \ , ° 2 9a . . 

-qkQi 
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From now on we will denote the relative acceleration of the contacting parti- 
cles at their contact points - the contact acceleration - as a a instead of g a . 
Equation ( |l6| ) turns into 

a a = b a + Aafsfp , (18) 

where again summation over (3 is implied. We will call this equation the ge- 
ometry equation. By means of this equation and the consistency conditions 
introduced in Sec. ^ we can now determine the contact forces fp. The con- 
sistency conditions read 

a a > 

U > (19) 
a a f a = . 

The first condition prevents deformation of particles, the second one excludes 
attractive forces and the third one requests that contact forces may only act 
if the particles stay in contact, i.e., if a a — 0. These conditions together with 
the geometry equation ( |T8] ) allows to determine the unknown contact forces 
f a . The whole system consisting of Eq. (|l8|) and the conditions ( |l9| ) is called 
a Linear Complementarity Problem. It can be solved by Dantzig's algorithm 

§■ 

To incorporate friction we introduce additional motion constraints which 
shall, if possible, impede a tangential motion of the contacting particles. For 
polygonal particles they are of the form 

g(q) = tj (n + Xi - rj - Xj) . (20) 

This constraint has a similar form as the constraint of the normal motion 
of the particles ([)]) but instead of the normal unit vector of the contacting 
face of particle j the tangential unit vector tj appears, thus ensuring that the 
edge of particle i at position r s ; + Xi does not move along the face of particle 
i away from the point rj + Xj of initial contact. 

These motion constraints are, however, of different nature than the normal 
motion constraints. Whereas in the case of the normal motion the constraints 
must never be violated, the constraints on the tangential motion may actually 
be violated, as it happens when the particles start to slide. This is due to the 
fact that the magnitude of the friction forces are limited by /i/jv , with fi being 
the friction constant and /n the corresponding normal force. As reflected by 
the consistency condition 4 (ref. Sec. |^) the friction force must adopt its max- 
imum value if particles actually slide, i.e. if the tangential motion constraint 
is inactive. Thus, in this case the value of the friction force is determined 
without need of further consideration. In contrast to the case of normal mo- 
tion constraints we may not neglect the inactive constraints because now the 
corresponding contact forces are non-zero and thus the constraint has to be 
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kept in consideration in order to determine the direction of the tangential 
force. 

Since friction causes only further motion constraints there is, in principle, 
no need of further discussion of the problem. The geometry equation ( |l8| ) can 
describe systems with friction as well. For simplicity of notation it is worth, 
however, to consider normal and friction forces and their corresponding mo- 
tion constraints separately. We now have 2 classes of motion constraints 



g N (q) >o 



, . ( 21 ) 



and the corresponding contact forces / and f T . Now the equation of motion 
reads 

For the second time derivative of the motion constraints we obtain 
■■n _d 9 % \ , dg» . ( v N dg» d$ 



(23) 



Renaming again ga N ' T ^ by aa N,T ^ and using the abbreviations 



oqk V > Oq k dqi 

oqk V > oq k dqi 



(24) 



and 



K ( M k Qk ) + tttsH** 



|«* 3 £«* ^ (25) 

^ " 9^ Mfc 9^ Aq/5 " Wk k Wk 
we can write the modified geometry equations 







(26) 
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The full set of consistency conditions then reads 



a. 




a 



■N 



(27) 



Equations ( |2q ) together with the conditions (|27]) can be solved with a modi- 
fied Dantzig's algorithm jl(| which will be discussed in the next section. Note 
that some of the tangential forces may be directly determined by the respec- 
tive normal forces, i.e., when sliding at this contact occurs. The consistency 
conditions for these forces have to be fulfilled, nevertheless. 

6 Dantzig's Algorithm 

Contacts can be classified into breaking contacts (a N = thus f N = f T = 
0), permanent static contacts, (a = and a T = 0) or permanent sliding 
contacts (a N = but a T ^ 0). If we knew a priory into which category 
each contact belongs the contact forces could be determined by solving an 
inhomogeneous system of linear equations which consists of all equations 
for which either = or = 0, with the corresponding f£ and /J 
as variables. All remaining normal forces are zero, the remaining tangential 
forces assume their maximum values. Unfortunately the contact classification 
is only known if we know the contact forces as well. 

We apply Dantzig's Algorithm to determine the forces together along with 
the corresponding contact classification. It starts with considering a certain 
contact, disregarding all others, i.e., their contact forces are set to zero. After 
having found a solution for this contact its classification is also known. The 
algorithm proceeds then with the next contact. Again the contact forces are 
determined preserving the consistency of all contacts considered before. In 
this process the contact classifications of the already consistent contacts may 
be changed if necessary. The process is repeated until the last contact has 
been classified. 

All contacts are assigned to one of four lists: 

• List NC of breaking contacts 

• List Cf of permanent static contacts 

• List of permanent sliding contacts. In list C + are all contacts where 
f T = fJ,f N , in C~ all contacts where f T — —fif N . 

All contacts in the above lists are considered to be consistent, i.e., they sat- 
isfy the consistency conditions. The classification is done successively for all 
contacts a by: 
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1. Check if the normal force = satisfies the consistency conditions 
a a > 0. If this is the case the contact is consistent and belongs to NC. 

2. If < we have to increase the normal contact force f£ to obtain a 
non-negative normal acceleration. However, increasing the normal force 
fa will change the contact accelerations of the already classified contacts 
as well. Since for the permanent contacts (3 we need dp = (and for the 
static contacts also aT = 0) this would invalidate the classification of 
these contacts. To preserve the consistency of the classified contacts we 
will have to change the contact forces of the the permanent contacts 
as well. We now have to determine how much we have to change these 
contact forces for a given increase s of the new force f£ to keep = 
and, if necessary aj = 0. To calculate the necessary changes we formulate 
a reduced set of geometry equations: 

n _ „N _ h N , A NN,rcd f N , A NT,red f T , A NN f N 

u - a f3 -Op + A pi j 1 + a 01 j 1 + Ji 0a j a 

n - n T -h T -4- A TN > Icd f N -4- 4 TT < rcd f T _|_ A™ f N 
U — dp —Op -t- J 7 -f- Jip 1 J 1 + Ap a J a . 

The reduced set of geometry equations can be obtained from the original 
geometry equations ( p6|) disregarding all breaking contacts and replacing 
fp — ifJ-fp for all contacts (3 from C ± (sliding contacts). Since only 
permanent contacts remain all contact accelerations in the reduced ge- 
ometry equations are thus zero. If we now change the new contact force 
f£f — ► f£ + s we have to vary the previously known contact forces fp N,T ^ 

in the reduced geometry equations by an unknown amount Afp ' in 
order to keep the contact accelerations at their required value of zero: 

=bf, + A™'™* (f» + Af») + A™ if? + Aff) 

\NN ( f N 



^ + A™' rcd (f, N + A f") + A 7, ,red ( A T + A fi) 

TN I rN , \ 



(29) 



Inserting Eq. (E8h we find 



=A™>™ d Af» + A^ d AfT + Afft 
=A T p^ cd A fj N + A T p^Afi + A™s 



(30) 



r jy rp\ 

This is linear system of equations for the unknown Af} ( ' with the 
step size s being a parameter. The necessary variations Af^ ' are 
proportional to the step size s, hence the solution is of the form 

Af$W = F^ T h , (31) 

where F-j is the necessary variation if s = 1 . By inserting the changed 
values back into the original geometry equations we now know as well by 
how much the accelerations of the breaking and sliding contacts change. 
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3. We increase the value of the new force f£? until either the acceleration 
a a = (the contact is now consistent) or until the classification of any 
already consistent contact changes. The classification changes if either 

(a) a normal acceleration a^j > becomes zero: The contact becomes per- 
manent and has to be moved from NC to either Cf or C ± according 
to its present value of the corresponding tangential acceleration. 

(b) a normal force fp > (permanent contact) becomes zero: the contact 
is now a breaking contact and has to be moved from Cj? or to 



(c) a tangential acceleration aT ^ becomes zero and the corresponding 
tangential velocity is zero: the contact now is static and has to be 
moved from C ± to Cf. 

(d) a friction force previously of smaller magnitude than its allowed max- 
imum reaches the maximum value ±/z/^ : the contact becomes sliding 
and has to be moved from Cf to C^. 

4. if we have not yet found consistent values for and we have to 
proceed with step 2. 

If the contact a is permanent we have to consider the tangential compo- 
nent of contact a too. The procedure is very similar to the calculation 
of ■ The only difference is that if > we have to decrease the friction 
force until it assumes its negative maximum value and if < we have to 
increase the friction until it adopts its positive maximum value. 

7 Collisions 

In the framework of Rigid Body Dynamics collisions occur if two contacting 
particles have a negative normal relative velocity at their contact point. We 
can easily convince ourselves that no finite contact force can prevent a de- 
formation of the particles. No matter how large the force is, it will always 
take a finite, however small, time to stop the approaching particles, hence 
they will deform each other. Thus, to prevent deformation of the particles an 
infinite repulsive force of infinitesimal duration is necessary. It turns out that 
the total momentum transfer Ap between the two particles is finitepf 



Resolving multi-particle collisions it turns out that the resulting state of the 
particles after the collision is not unique. This is due to the infinite stiffness 

2 Note that the duration t c of the collision is only formally the parameter of the 
limit in the equation above. This limit is actually achieved by starting with 
deformable particles and increasing their stiffness to infinity. In this limiting 
process the duration of collision approaches zero. 



NC. 




(32) 
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of the particles, or equivalently, the infinite speed of sound in the particle 
material. In the limit of infinite speed of sound all information on the exact 
collision mechanism, e.g., the sequence of individual pair collisions, is lost, 
since any of the collisions is of vanishing duration. Colliding particles of finite 
stiffness do not exhibit this feature, since here all processes take finite time. 

According to these arguments the necessary information on the detailed 
collision mechanism is not accessible. The following set of assumptions turned 
out to yield realistic results, although they cannot be uniquely derived: 

1. All individual pair collisions occur at once. 

2. The transfer of momentum at contact points is finite. There is no mo- 
mentum transfer which corresponds to attractive forces. 

3. The relative velocity after a collision can never be smaller than — ev, 
where e is the coefficient of restitution and v is the impact velocity at 
that contact (v < !). 

4. If the velocity after the collision is strictly larger than — ev there is no 
momentum transfer at this contact. 

The first two assumptions have been introduced already. The remaining two 
assumptions deserve further discussion. Two-particle collisions can be de- 
scribed by means of the coefficient of restitution which relates the precolli- 
sional relative velocity v and the final velocity v' after the collision: 



In the case of multi-particle collisions, however, two particles which initially 
rest relatively to each other may separate after a collision. Therefore, the 
final velocity may indeed be larger than the value — ev. Contrary it may no 
happen that two particle which collide with a finite impact velocity are at 
rest relative to each other afterwards. Therefore, v 1 > — ev. 

The fourth assumption simply means, that if two particles separate from 
each other with higher velocity than — ev their aftercollisional velocity may 
not be increased further by an additional momentum transfer. 

For convenience we define the excess velocity Av 



which is zero if v' = —ev. The source of any velocity change is a momentum 
transfer between the colliding particles. We can relate the excess velocities 
at the contacts with the momentum transfers by means of the collisional 
geometry equation, which is derived in a similar way as the geometry equation 
of the force algorithm: 



v' = —ev . 



(33) 



Av = v' + ev 



(34) 



Av a = (l + e)v a + A$A Pf) . 



(35) 



Note that the v a are relative velocities of the contact points of two contacting 
particles, but not the velocities of the particles themselves. In mathematical 
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terms the above discussed assumptions read 

Av a > 

A Pa > (36) 
Av a Ap a = . 

The first condition prevents the final velocity from being smaller than — ev. 
The second condition excludes attractive interaction between particles. Fi- 
nally, the third condition means that there may be a finite momentum transfer 
only if Av = 0. These conditions together with the geometry equation ( ^5|) 
form a Linear Complementarity Problem which is already familiar from Sec. 
H and can be solved by Dantzig's Algorithm. 



8 Resolution of static indeterminacy 

If the number of contacts in the system is too large, the contact forces cannot 
be uniquely determined by the force algorithm. This occurs, if the number of 
free variables in the system, i.e. the number of contact forces, is larger than 
the number of mechanical degrees of freedom, 3N in 2d or 6N in 3d, with N 
being the number of particles. However, the total forces and torques acting on 
the particles and, hence, their trajectories are unique. This drawback restricts 
the applicability of Rigid Body Dynamics for the simulation of railway ballast, 
since this system is of mainly static nature, the exact knowledge of the contact 
forces is crucial for understanding its behavior. 

So far we have considered the contact forces as independent of each other. 
This assumption is the reason for the force indeterminacy. In realistic systems, 
however, the forces are not independent (Fig. ||). If we apply an external force 
directed to the right on the central particle, we increase the contact force with 
the particle to its right while at the same time decreasing the contact force 
with the particle to its left. In this example both contact forces in reality 
depend on a single parameter, which is the applied external force. 




Fig. 3. The central particle is 
forces on both contacts are not 



in contact with two other particles. The contact 
independent of each other. 
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We can mimic this behavior by introducing small displacements of the 
particles. Each particle has a set of macroscopic coordinates r and <p an d a 
set of microscopic coordinates Sr and 5<p. In analogy to the vector q of all 
(macroscopic) coordinates of all particles, we define the vector Sq of all mi- 
croscopic (infinitesimal) coordinates. The contact network and the kinematic 
state of the particles is determined solely by the macroscopic coordinates, 
hence they can be considered as the actual coordinates of the particles. The 
microscopic coordinates lead to a deformation of contacting particles. By def- 
inition they are of infinitesimal size, which allows us to restrict ourselves to 
a linear approximation in 5q for the computation of the deformations. The 
vector £ of all deformations in the system 

£ = DSq (37) 

is defined by the deformation matrix D and the microscopic coordinates. The 
dependence of D on the geometric properties of the systems is straightforward 
but lengthy, therefore we will not give explicit expressions here (for the full 
derivation see fl5fl )- We define a force law to relate the deformations at the 
contact points with the contact forces: 

/ = f(0 ■ (38) 

Now the contact forces are functions of the displacements Sq. To calculate 
the forces we use the geometry equation (^6j) together with the consistency 
conditions (p7|). This set of equations is to be solved for the microscopic 
displacements 5q. Hence, there are as many variables as degrees of freedom, 
i.e., the system is unique. 

To determine the microscopic coordinates we use an overdamped relax- 
ation method. We start with a set of inconsistent coordinates 5q. Inconsistent 
means that the consistency conditions for the resulting forces and accelera- 
tions are not fulfilled. Now we let the system relax. This way we find new 
microscopic coordinates such that £' = £ — ha, with h being the step size 
and a the vector of all contact accelerations. If the contact acceleration is 
negative the displacement will be larger, yielding a larger force to stop the 
approaching motion of the contacting particles. Since the adjustment step for 
the displacements is proportional to the acceleration the sequence of micro- 
scopic coordinates of the particles can be understood as overdamped motion. 

The adjusted microscopic coordinates are the solution of the linear system 
of equations 

£-ha = DSq' (39) 

or, equivalently 



D (Sq' - Sq) = -ha 



(40) 
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In cases where there are less degrees of freedom than contact accelerations 
the system of equations is overdetermined. There may be vectors ha which 
are not representable by any vector Sq' — Sq. In this case we have to project 
the vector —ha into the image space of the operator D before solving the 
system of equations. 

We repeat the adjustment of microscopic coordinates until the consistency 
conditions are met. To further improve the speed of this method we can save 
the microscopic coordinates which yielded a consistent system in the previous 
time step. If the system did not change too much this set of microscopic 
coordinates is very close to the new solution and we need only few iteration 
steps to arrive at the new solution. 

This method combines the advantage of Rigid Body Dynamics, namely 
the ability to simulate very stiff particles, with the advantage of Molecular 
Dynamics, namely uniquely defined contact forces. An additional advantage 
of the method of small displacements is that we can now easily simulate cer- 
tain degradation mechanisms. If we, for example, want to simulate the effects 
of abrasion of edges of the particles we can do this by gradually changing the 
force law (p8|), which describes changing edge properties of the particles. 



9 Step size control 



The integration scheme used in our implementation is a Runge-Kutta method 
of fourth order. During one time step there are four force computations nec- 
essary. Since we use discrete time steps we are frequently faced with the 
problem that after a given time step some particles do deform each other. In 
this case, obviously, the chosen time step was too large. This problem could 
be solved by predicting the time of next contact (collision) from the present 
state of the particles. Since the particles are subject to forces which vary in 
time this prediction cannot be accurate, as we approach a collision we would 
have to update it repeatedly. Since this prediction method is quite compli- 
cated we chose a simpler method. We advance the system by a certain time 
step. When determining the contacts in the next time step (see Sec. ||) we 
check for deformations of the particles. If any deformations occur we restore 
the state of the system before this time step and take a time step of half its 
original value. We then repeat this process, i.e. advance the system by the 
new time step and check for deformations. If we again encounter deformations 
we again divide the time step by two and repeat the computation until a state 
without deformations is reached. The new state of the system is accepted and 
the other steps of the algorithm are performed. This procedure ensures that 
there will be no particle deformations in the system at the beginning of any 
accepted time step. 
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Since this method can only decrease the time step, however, we need a 
procedure to increase the time step again^] in order to avoid permanently 
slowing down the simulation. We can always increase the time step if there 
are no collisions in the systems for a certain time. Hence, if we successfully 
performed a number of time steps without encountering collisions we can in- 
crease the time step again by a factor of two. We cannot increase the time 
step after only one successful time step since this may lead to frequently al- 
ternating increase and decrease steps. A requirement of three successful time 
steps before increasing the time step has shown to yield good computational 
efficiency. 

With this choice of a step size control method we have completed the 
discussion of the simulation algorithm. 

10 Conclusions 

Although Molecular Dynamics methods have recently been applied success- 
fully to the simulation of the dynamics of railway ballast [jl6| we have doubts 
that this technique is suitable to the simulation of almost rigid, sharply edged 
particles such as railway ballast, in particular if the system behavior is gov- 
erned mainly by static properties of the system. Moreover, it seems to be 
unsuitable for the simulation of long time effects such as densification and 
wear of ballast. 

The Rigid Body Dynamics is a method which is intended to describe the 
motion of systems of very stiff particles. In a natural way it avoids the prob- 
lems of Molecular Dynamics simulations which have been discussed in detail 
in Sec. ||. These problems are unavoidable within the concept of Molecular 
Dynamics. Therefore we believe that Rigid Body Dynamics is much better 
suited to the simulation of railway ballast than Molecular Dynamics. From 
the presentation of the algorithm we have seen that a time step in Rigid 
Body Dynamics is a lot more complicated than a time step in Molecular Dy- 
namics. Hence the according implementation requires by far more computing 
time than an implementation of a Molecular Dynamics code. Fortunately this 
disadvantage is offset by the fact that in Rigid Body Dynamics we can chose 
a larger time step than in Molecular Dynamics. Whereas in Molecular Dy- 
namics the time step is determined by the critical deformation (the length on 
which the force changes by one unit) divided by the characteristic velocity in 
the system, which usually yields a very small time step (e.g. ~ 10 -7 sec), the 
time step in Rigid Body Dynamics is determined by the characteristic time 
in which the forces in the system change significantly. The latter quantity is 
usually much larger (~ 10~ 3 sec). Thus, the computing time requirements 
for both methods are comparable when simulating mainly static problems as 
it is the case for railway ballast. 

3 There is, of course, an upper limit to the time step which is dictated by the 
accuracy of the chosen integration scheme. 
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An additional advantage of Rigid Body Dynamics is that fracture of par- 
ticle can be introduced in a very natural way. When dealing with dcformablc 
particles in MD simulations the fracture of particles can lead to numerical dif- 
ficulties. After the fracture of a particle, i.e., when a large particle is replaced 
by two or more fragments, there will be a finite time in which the system 
has to readjust and to reach a new stable state. There are situations when 
the new stable state deviates significantly from the realistic (experimental) 
situation. This problem is not encountered in Rigid Body Dynamics. After 
a fragmentation (which is algorithmically done by constructing two or more 
new polygons from the previous particle simply by introducing another edge 
into the old particle and splitting along this line) the system remains always 
stable. 

It is our firm believe that Rigid Body Dynamics provides a very promising 
alternative to Molecular Dynamics for the simulation of railway ballast. 
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